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A  mJMERICAL  METHOD  FOR  THE  RADIATIOM  TRANSPORT  EQUATION 
IN  PLANE  GECMETRY. 


I.  INTRODUCTION 

The  scope  of  work  reported  here  is  the  development  of  a,  numerical 
■  method  for  performing  time  dependent  calculations  of  radiative  energy 
transport  through  slahs^  Or  guiding  interest  has  teen  the  retaining  of  • 
contact  with  the  diffusion  approximation  hy  maintaining  a  parallel  of 
numerical  methods.  In  doing  so,  it  has  "been  found  that  the  radiation  ' 
intensity  may  effectively  he  eliminated  from  the  equations  and  from 
computer  storage.  As  is  the  case  in  the  diffusion  approximation,  only  the 
material  temperature  remains.  Conservation  of  energy  is  also  a  direct 
consequence  of  the  method.  However,  to  accomplish  this  objective  a  number 
of  simplifications  have  been  made  which  may  limit  the  possibility'  of 
generalizing  the  method.  These  include  the  assumptions  of: 

1.  Plane  slab  geOT(etry  and  boundary  conditions  of  several  simple 

types.  y. 

a.  Free  surface  radiating  into  vacu\mi. 

b.  Opaque  surface  with  zero  flux. 

c.  Specified  forward  or  backward  flux  into  the  slab. 

2 .  Negligibility  of  the  time  for  photon  flight  between  emission 
and  absorption.  Related  neglect  of  the  specific  heat  of  the 
radiation  field  in  comparison  with  the  specific  heat  of  the 
material . 

3*  Local  thermodynamic  equilibrlian  of  the  material.  ■ 

4.  Interaction  of  radiation  with  material  only  through  true 
absorption.  No  photon  scattering.  Frequency  independent 
absorption  coefficient. 
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Several  of  these  assumptions,  namely  2,  3,  and  neglect  of  scattering, 

■  may  be  satisfied  by  restricting  consideration  to  low  temperature  transport 
through  materials  of  not  too  small  density.  The  more  serious  approximation 
of  frequency  independence  of  the  absorption  coefficient  has  been  made  in 
order  to  separate  the  two  important  problems  of  the  geometrical  behavior 
and  the  frequency  behavior  of  the  radiation  transport  equation.  .  Extension 
of  the  present  method  to  .include  groups  of  photons  having  different  absorption 
coefficients  will  probably  not  offer  serious  difficulties-  to  the  present 
method.  However,-  it  is  by  no  means  clear  that  extension  of  this  basic 
procedure  to  other  geometries  can '.be  made.  In'  this  sense,  the  present  work 
is. of  interest  in  its  special- applicability,  to  slab  problems.  In  compensation 
it  offers  a  rather  simple  system  of  equations  having  small  storage  require¬ 
ments-  and  light  stability  restrictions  which  allow,  rapid  solution  of  radiation 
transport  problems.  ' 

■  ■  II.  IWTBGRAL  EQUATION  FOR '  RADIATIVE  TRAKSITO  . 

-.  •  The  assumptions  of  plane  geometry  and  azimuthal  symmetry  of  the 
radiation-  intensity  permit  considerable  simplification,  of  the  equations  of 
radiative  transfer  i  As  will  be  shown,  the  equation  for  the  rate  of  change 
of  material  ener^  for  this  case  can  be  reduced  to  an  integral  equation  in 
which  all  reference  to  the  intensity  of  the  radiation  field  has  been  elimi¬ 
nated. 

The  equation  for  the  rate  of  change  of  the  intensity  of  the  radiation 
field  I-i;  is: 

31^  di-o'  ■  ■  . 

aST  "  ^  -  .  .  (1) 

4  =  cos  ©  ,  ©'=  polar  angle  of  radiation. 


where 
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■  '  hj>  ' 

.  .  =  (l-e  )  a~ )  cr~  ~  ateorption  coefficient  and 

<7:5'  takes  account  of-  the  induced  emission  of  the  material.  . 

_  2h 

-  1 

is  the  Planck  function  source  assigning  that  the  material  collision  rate  is 
sufficient  to  establish  local  .thermodynamic  equilibrium  in  the  material 
cheuracteriz'ed  .by  a  material  temperature  T .  Photon  scattering  and  polari¬ 
zation  are  neglected  in' Equation  (l). 

Dropping  the  subscript  indication  of  frequency  and  introducing  the 
optical  depth 


the  intensity  I  at  position  -jr-and-in  direction  ^  can  be  expressed  from 
Equation  (l)  in-  terms  of  the  Planck  function  at  points  along  the  photon 
pencil  of  radiation: 


rr  .r  S’ 

B(?r^)  e»"  d-C*  +  I  (r;) 


‘-o 

The  effect  of  retardation  or  time  of  photon  flight  from  the  so\n:ce  has  been 
neglected  so  that  the  intensity  is  approximated  as  depending  only  on 
instantaneous  values  of  the  Planck  function.  This  corresponds  to  neglecting 
the  first  term  of  Equation  (l)  and  is- justified  -when. the  temperat-ufe  -waves 
travel  much  slower  than  c.  The  integral  extends  into  the  material  to  an  • 
arbitrary  depth  at  which  point  .the  intensity -in  the  n  -  direction  is 
I  This  point  may.  be  taken  as  the  boundary  of  the  system  or  it  might 

be  a  depth  beyond  which  the  contribution  to  I  is  negligible . 


Having  neglected  the  time  dependence  of  the  Boltzmann  equation,  the 
dependence  on  time  oif  the  radiation  transport  equation  arises  from  the  change 
of  energy  as  governed  by  the  first  law  of  thermodynamics.  In  general,  the  • 
total  Internal  energy  changes  as  a  result  of  motion  of  the  material  and 
of  entropy  changing  processes  such  as  heat  conduction.  We  wish  to  set  in 
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evidence  In  this  discussion  the  characteristics  of  radiative  transfer; 
consequently,  we  neglect  all  other  mechanisms  of  energy  change,  including 
work,  and  write: 


5i-  ^  ^ 

at  ax  at 


(3) 


where  is  the  material  internal  energy  per  unit  volume,  is  the  radiant 
energy  per  unit  volume  given  by: 


;  =  2* 


'OO 

dV 


I  <4i 


and 


F  =  2jrc 


c  f  I 

^  -^-1 


I  jid+i 


.(J+) 


is  the  component  of  radiant  energy  fltuc  along  the-x-  direction.  Since  a-', 
numerical  approach  paralleling  the  diffusion  treatment  is  ultimately  desired 
.Equation  (3)  is, integrated  over  a  small  Interval  to  give  the  Integral  ’ 
expression  for  energy,  conservation.  Indeed,  the  integral  equation  could 
equally  well  have  formed  the  starting  point  of  the  investigation.  In  the 
res\iltlng  conservative  form  the  way  to  insure  conservation  of  energy  in  the 
difference  equation  becomes  cleeur.  The  desirability  of  this  conservative 
property  of  the  equatipns  in  fact  dictates  the  choice  of  the  following 
system  over' alternative  .forms  based  on  the  differential  equation,  Equation 
'(3).  ■■  ■  ■ 

.Xg  .  _  • 

(5) 


f  dx  =  F(x^)  -  FCxg) 


jC 

■^1 

in  the  above  we  have  neglected  the  rate  of  change  of  radiative  energy  in 
comparison  with  that  of  the  material  energy  as  is  Justified  for  lew- 
temperature  Interactions . 

In  order  to  remove  the  explicit  appearance  of  I  from  Eqmtion  (4) 
for  the  fluxes,  we  use  Equation  (2)  to  obtain  an  integral  expression  in 
which  the  sources  of  radiation  in  the  emission  of  the  material  are  displayed 
Forward  and  backward  directions  are  treated  separately:  for  the  forward 
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direction,  |i  =>  0,  the  integral  extends  from  =  -oo  where  I  =  0; 

for  \i  ^  0  we  take  =  +ooand  again  assume  I  (z:^)  =  0* 


'  -  Zzl!. 

f 

I  (idji  =.|-  B  dr*  I 

e  ,  ^  •  dji 

Jq 

-ico  ^0 

r°  ■ 

rOO  • 

.  zzzL 

I  ^ldtl  =  B  .dZ  r 

(-l)e  K 

u 

■ 

(4i  .  (6) 

•  The  angular  integrations  in  Equations  (6)  are  particularly  simple 

because  the  source  B  does  not  depend  on  the  angle.  Introducing  the  E  (x) 

2  ■  ■  ■■  '  ■  ■  •  • 
functions  defined  by:  •  •'  •  =  .  •. 

X.  •  •  •  '.v 

E  (x)  =1  e  d|i  .  ,  '  • 


these  expressions  become:  • 


. .  ■  ,1 


'Zr  . 


,  I:  ^d^  =  ' .  bCc'  X (Tr •  y-dV  ■  ■■,  ■■ 


^OO 


Equation  (5)  becomes; 
•  r^2  - 


rco  ■ 


dx  =  2nc  ■■■' '  d^ 

:i  ..  .  ■  "■  -Jp-  ■  . 

+  [..  'b.>  E^.- (V  -  dV 


r^i  •  ■ 

+  B 


.  •  •  ^co. 

B  vie^  (?r^  - z:')  dT^  -  .  b  •  (z^.  -  q  )  dz:' 


(7) 


■-  .  B  '.q  (q, :  z")  dZ  ..  ■ 


(8) 


Equation  (8)  which  is  the  starting  point  of  the  numerical  work  expresses 
the  dependence  of  the  rate  of  change  of  energy  in  the  selected  region  on 
the  attenuated  sources  as  they  compose  the  radiative  fluxes.  The  appearance 
of  functions  instead  of  the  exponential  arises  from  weighting  with  more 
■5 - 

see,  for  example,  Appendix  A  of  "Introduction  to  the  Theory  of  Neutron  . 

Diffusion"  by  K.  M.  Case,  George  Placzek,  F.  deHofftnann,  Los  Alamos 
Scientific  Laboratory  notes. 
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strongly  attenuated  otligue  rays  while  the  positive  and  negative  terms 
correspond  with  forward. and  backward  contributions  to  the  flux  Integral. 

We  illustrate  the  nature  of  the  contributions  to  Equation  (8)  by  schematically 
indicating . the  behavior  of  the  kernels  of  the ■ integrands .  ■  ’ 


For  completeness,  we  return  to  the. differential  expression, of 
Equation  (3)  to  tabulate  several  alternative  forms  which  may  be  of  interest. 
Forming  ^  from  Equation 


in  which  Equation  (l)  has  been  used.'  Equation  (9).  expresses  the  statement 
that  the  rate. of  change  of  material  energy  is- the  difference  between  rates 
of  absorption  of  I  . (which  occurs  only  in  the  average  Intensity)  Equation  (2) 


for  I  is  integrated  separately  in  the  forward  arid  backward  directioiis. 
Thus,  the  integrals  contributing  to  E^,  similar  to  Equations  (6),. arris 


('■ 

1 

^  ■ 

r'C 

I  du 

=  f  B 

•  dr’/ 

e 

d|i 

=  f.  . ,B.-.  •. 

■  J-oo 

Jo 

|i 

J-OO 

.00 

A 

rOO  '. 

I  qi 

.  .B- 

■  -"z: 

.d^l- 

J'o 

e 

:  .=' 

;'.B  >  E 

'ZJ  ■- 

Using  these  results  in  Equation  (9),  the;  equation  determinirig  the.  temperature 
of  the  material  is;  •  ' 


g  =  2jtc  .  C7^  dl/l  +  ..  bC?:’  )  E^  (ir  -  I)  dr*- 

■  Jo  L  .  J^oo  ■ 


(10) 


where  the  specific  heat  at  constant  volume  of  .  the  material  and  B.('2r) 

depends  on  position  through  the.  temperature.  •  Use.  of  Equation  (2)  for  I  in 
which  photon  retardation. has. been  neglerited  implies  in  Equation  (lO)';'  that 
the  energy  in  the  radiation  .field  must  be  small'  eno\^h  compeured  with  that  of 
•  the  material  to  warrant  omission. .from  the  energy  equatiori.  .  Equation  (10 ) 
is  a  nonlinear  integro- differential  equation  for  the  temperature  containing 
no  further  reference '.to.  the  radiation.:  field  in'tensitles.  The  terms  correspond¬ 
ing  to  the  emission  and- absorption  of  energy  by  the  material  are  displayed 
as  the  negative  and  positive'members  of  Equation ■■( ID).  The  limiting  case 
of  optically  thin  material  in  a  vacuum- is  given  by-  the. first  term;  the 
frequency  Integral  of  it  gives  the  well  known  "Planck  Mean  Absorption 
Coefficient"  .'which  is '  frequently -cbraputed  .-together  -with' the  ^Rosseland  mean  '  ..  ... 
from  absorption  coefficient  -values.'  ■  ■ . 

The  integral  describing  the  cbntribution  from  absorption  tP  ithe.  •.  .-..  -? 

energy  change  in  Equation  (lO)  can  be  transformed  thr&ugh' .integration  by  ^ . 
parts..  '  ...-  ^  .  i' ,  ■  ■  ’ 

f  B(r')  E,  (ir-  D'dr  =  2  B(r)  (r'  -  ^  ^ 

jLoo  ... 
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Substitution  in  Equation  ( 10 )  gives  the  following  alternative  form: 


2jtc J  (T*  dt)  I^J  g,  Eg  (r'  -  r)  dr 

H  .  •  -v 


''v  dt 

'Zr 


Further  integration  by  parts  in  Equation  (ll)  yields: 

■  „  ■■  f'Oo  ,  ,00 

•  .'C  '  H  =  2itc I  dV  -±±- E  or  -  V  i)  dr* 

•  7  Jo  •  .  j-ood(-rr  ^ 


(ii) 


■  (.12.) 


Ill .  LIMITING  EXPRESSIONS  ANIl  APPROXIMATIONS'  .  ’ 

Equation  (12)  is  a  convenient  form  to -discuss  the  limiting  case-  in' 
which  the  temperature  changes  by  a  small  amoxint  within  a  mean- free  path. 
Expansion  of  B"'  (r' )  about' the  point  -er  allows  the  evaluation  of.  the  lead¬ 
ing  terra  of  the  Integral  in  Equation  (12).  The  result  is:  ’ 


1  ^ 

"v  dt 


lute  kne  d  r'^dB  dV 

~  H  ^  ^ 


where 


,00  ■  •  L 

4«c  d  /dT  I  A.  f  ^  ) 

3  dx.  W  Jo  ■  dT  '  3  dx  dx-  ' 


I  ^  dlJ:' 

Jo  dT  ^ 


R 


.  ''(IS) 


/•  ^  dll  •  ■  ■  .  .  .  ■ 

is  the  Rosseland  mean- free  path.  Equation  (13)  gives  the  radiative  diffusion 
approximation  in  which  the  energy  of  the  radiation  field  has  been  neglected. 

The  optically  thin  result  obtained  from  Equation  (lO)  is  the  opposite 
extreme  in  approximation 
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dT  4 

C  -  (T  ca  T 

V  dt  p 


(li^) 


where 


f 


B  cr^’  6n) 


(^■p  ■  ■■  ^co 


■B'd'l^ 


Je¬ 


ts  the  Planck  mean  absorption  coefficient. 

If  the  absorption  coefficient  depends  oh  type  of  materiaiy' density 
and  teini)erature  bat  not  on  the '  frequency  (the '-."gray  atmosphere"  ■appr9ximation). 
the  integration  over  frequency  can  be  carried  out  -in  Equation  (lO);  '  , 


-(I'T' 

■C  -Tz  -  CT"'  ■  ca 
V  dt 


1  I  fn^ 


-  (Z)  *  (V  )  -  Z  1). it! 


■(15) 


The  corresponding  approximation  to  Equation  (S)  is  obtained  by  using  the 
following  value  of  F(x^)  in  Equation  { 5 ) :  ' 


F(x-  )  = 


ca 


■y  2 


J  'ey  -iV  T^(t:' )  .1^  (r'  -■fj )  i't 


.(16.) 


IV.  DIFFERENCE  EQUATIONS ’ FOR  RADIATION  TRANSPORT 


In  formulating  a  suitable  numerical  approximation' to  Equation  (5) 
we  keep  in  mind  several  criteria.  The  method  should 

■1.  be  as  simple  and  easy  to  compute  as  possible 

2.  retain  property  of  energy  conservation 

3.  be  sufficiently  acevu-ate  to  limit  to  the  diffusion  and 
transparent  regimes  properly 

4.  have  an  adequately  wide  region  of  numerical  stability  to 
permit  fast  integration  forward  in  time. 


In  conformity  with  the  radiation  diffusion  -  hydrodynamics  treatment  j, 
the  region  of  interest  is  divided  into  a  limited  number  of  adjacent  slabs 
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or  zones .  In  each  zone  the  characteristic  temperature  is  chosen  so  that  . 
the  corresponding  specific  energy  vhen  multiplied  hy  the  zone  mass  gives 
the  correct  total^  zone  energy.  As  a  first  approach  to  this  method,  we 
devise  a  simple  explicit.  ■' treatment  of  the,  temperature.  The  time  derivative 
of .  the  material  energy  is  written  as:-  . 


jVl/2.  VJ+1/2  >^n+l  _ ,  n  , 

■  .  •At.;',  ^  .V  j+l/2  ■  J+l/sA 


(IT) 


an  express,i.on  which  permits -  the  ;,temperature  .to  .he  explicitly  determined 
when  the  right  ■'hand  ,, side  of  ';  the  equation  :is  evaluated' with'- quantities  known 
'  at- ’time' n.  '-The'  q'uantitles-:,m^,|,'^^2,,C^j-Tj^^2  •,are  the''  zone  mass  and  -specific 
heat;  lat'tei;';'4uhnt'ity;''may.i.'depehd;on. time 'through  its' temperature  .  •  '  '  ’ 

■,'dependence>.’.,,,'As':'ihSl'cated'''-lh  a,', shhsequent- section,  'this' formulation  Is  ' 

,  e,xjpect,e:d-;;tp;,,'haye;;C^  stability.  Thefe  also  will' be  errors  of  second- 

order,'.' ini,  the’ iime'iihW  •,case''.in'’the'  corresponding  explicit  treat- 

ihent,vo,f-i-^he'|radia'tion;.,dl,ffu’si^  •  . 

-  r;In  -,cohfq:nmtt’'y^  above. -it  .'is  necessaiy  to  use 

some  careiihVhyaluati^,.!^,.^,'!!!^^]^^^!’? '.’contributing-  to: the  radiation  flux. 

In  ’par'tiGular, ;  ihe,:simple,,s'tr,s  in  which  the  source-  is  .taken  to  be  constant 

in  each  zone ’’does  hot  limit' properly  -to  the  diffusion  expres.sion.  As  would  • 
be  expected  from  .the-.-fact  ’tha't’' the  diffusion  flux'  is. given  as  a  first 
derivative, '-'it  is  necessary  in  the  -zones  immediately  adjacent  to  the  boundary 
at  which  the  flux  is  evaluated  to  use. a  "linear  interpolation  of  .the ; source 
'..strength.  In' zones  farther ’ away  the  criterion. of .  simplicity  su^'ests - 
assvunption  of  constant  source  strength  .for -the  "tail"  of  the-  integration. 
Since'-.the  independen't  variable  in.  the  flux  integrals'- is .  the  optical  depth  .  '  • 
-if  is -necessary  to  evaluate  it  at  zone  interfaces 'on  each  cycle 'of  the.  . 
calculation.- '  At  .time  n  the  depth -is  given  siiiiply'by: 

'•-  j-1  •'  ■'  . 


'■  -•■'T. 


'J' 


'd-' 


i+i/2  ^^i+l/2 


(18) 


in  which  cr’  is  evaluated  from  the  zone  temperatures  optical 

depth  of  the  zone  center-  ~  (^j  ^"^+1^’  Intervals 
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in  the  optical  depth  In  general  are  expected  to  d.epend  on  position  requiring 
that  the  interpolation  formiila  be  correspondingly  more  general.  • 

In  evaluating  F(x.)  from  Equation  (l6),  the  use  of  linear  interpolation 
is  restricted  to  the  half-zone  on  either  side  of  the  zone  interface.  In 
this  region  the  fourth  power  of  the  temperature  is  given  by  interpolation: ,, 


where 


■  (19) 


P 


/  4  ■  '  4 


'^d-1 


.1-1/2 


), 


^J-1 


Carrying  out  the  integration  introduces  E^.'for  the  constant  regions  and 
both  E^  and  Ej^  for  the  linear  regions.  The  first  integral  corresponding 
to  the -forward  current  in  Equation  ’(l6)  is.  given  by  the  sum  of  three  terms; 


l"j  =  “J 


■1/3 ■ 


/^:-'^-i/2^  ^3  ^ -:^^i/2^ 


2-j  “  -3-1/2 


^3  "  ^-1/2^  "  ^3  ^^3 


(20) 


3-2 

3^3  "IJ!"  '^i+1/2 

-oc 


corresponding  ■  to  the  linear  region /from  1/2'  constant  temperature 

Half- zone  .from  to -and  the  remaining  zones' of  constant  temperature 
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The  analogous  set  of  three  terms  for  the  hackward  curTent 

"  “j  ^”^+1/2  “ 

-  Pj  [V3  -  -  2^)  -.  (2rj,,/2  -  rj)  E3  -r^)]  •  .  , 

-"  -j+i/^  -  .®3  -(’^+l/2  ■  2^^  ■  ®3 
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.  are  to  he  summed  and  subtracted  from-.’tlie  forward  current  to  form  the  net'  . 

flujc  at  .'the  point  i  -.The-  complete  difference  equation' for  the .  temperature  ■ 
■  results'  .  '  '•  ■  ■.  ..  '  "■•■•  •  .■■  ■  •  ■  '  .  '• 
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.;  ;,from-  substituting  such  expressions  for  the  z one "boiindary 'fluxes  -in  the' energy  ■ 
Equation  ■( 5 )  of  which  Equation  ( 17 ). is  the- approximation  for  the  left  hand 
side'.  This  ■  expression  may  be  solved  for  i- '  ’  Ail  of’ the  contributing 

terms  are  assumed  to- be  function's -'of-- the  known  temperatures  ■ 

•may  be  replaced  in  storage  by  .the  .newly  calculated -values  at  the, ’.advanced 
--  .time.  Consequently,  the  method  may  be  termed  "an  explicit  integral  formu- 
•  -lation  of  the'  time- dependent  radiation  transport  equation". 

'  .  ’  Boundary  conditions  .-are  derived  by'  ,  considering  the  flux,  expression. 

A  free  surface,  condition. corresponds  to  setting  the  source  to  zero  'for  all 
zones  beyond  the  free  interface.  A. '"hot  wall"  bo\mdary  to  simulate  the 
problem  of  the  development  of  the  diffusion  wave  from  a  hot  source  corresponds 
to  giving  all  zones -beyond  a  given  .boundary  the  same  temperature. 
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STABILITY  OF  THE  DIFFBRSMCE  EQUATIONS 


The  system  of  difference  equations,  Equation  (5),  with  boundaiy 
conditions  and  initial  conditions  may  be  solved  for  the  temperatm'es  of  the . 
material  in  an  explicit  step-by-step  numerical  calculation  to'obtain. the- 
■solution  at  any  desired  time.'  In  addition  to- the  question  of  the  truncation 
error  incurred  in  the  numerical  approximation  the  question'of  stability  of 
the  solution  must  also  be  raised.  .  If  small’' errors  are -introduced,'  perhaps, 
througli ‘rounding  of  the  numbers,  if  is  necessary  that  the- difference  equations 
■not'  amplify  the  err.or.s  without  bound.-  It  is’frequent-ly  ’the' case  that  an, 
inequality  must  be  sa-tisfied  rest’iricting  the  size  of -the- time'  step  to  Insure 
the  stability  of  the  solution’ in -this  sense.’’  .  ' 

Investigation  of  stabi'lity  .of  the,  nonlinear -'equa-bion  with  boundary;.'  - 
conditions' is  a  job' beyond,  the  scope,  of  , this  Investigation.  ■  However,-'a',’.’ " 
'more  modest  effort 'will ’  indicate  some 'df  the -practi’cal  limitations  to.be 
■observed  in- the  numeric.al’ work  to  follow.  ’’As  is -generally  done  in'  this',  . 
situation,  the  equations  are  .linearized  and' the 'Von  IJeumann  criterion'is’  '. 
invoked’.  Introducing  the  .depende'nt  variable  fj+'2_’y2  “  the  equat’ipns -. 

become  linear  with  constant  coefficients  provided  that-, the  solution  is '■ 
restricted  to  small  variations.  ■  ■  '  ’’  "  " 


At 


'^^J+1/2  ^j+i/2”l  ’  ■  '2  ■  ■^'’,2^j  ■  ,3^j 


■  l^j+1  ■  2^j+l  "  ■  3^j  ^■•l'^j+1.'^ ’2.^jil. 


(23) 

The  f's  and  b's  are  functions  of  the  accordance  with  the  explicit 

forniulation  of  the  system.  Since  these  equations  are  now  linear  equations  ,, 
with  (asnujiied  locally)  constant  coefficients,  they  also  have  the  same  form 
as  the  equations  governing  propagation  of  small  errors.  The  solution  is  • 
assumed  to  have  the  form: 


J  (r,t)  = 


ik'T'  at 
e  e 


(24) 
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in  which  the  wave  number  k  plays  the  role  of  a  free  parameter  while  a  is 
to  be  determined.  For  stability 

(^aAt|  ^  ^  k  , 

since  the  above-  is  the  Mplification  factor  of  the  solution’  in  one-time 
interval . ■ 

•  In  order  to  further  simplify  the  analysis  the  rnesh  interval  and  . 
absorption ‘coefficient  are  assximed  constant.  Thus,  ■ 
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l^J  "  “j 


1/2  -  I!,  (§)] 


"^0 


1/3  -  Ej^  (f )  -  f  E^  (f ) 


2^J  =  (g)  -  E^  (a) 


.  oo 


iknA' 


n-x 


(nA)  -  E^  (n+1)  A 


in  which 


aj  -  1/2- (1  + 
■  (1 


Also 


_  -1^^  -fi  .  .  f».  -  f '  •  •  f'  '  =  '^f '  ■ 

rj+1".®  l^J  ’  2^J+l  “®  2^J  '  -3  J+l-  ^  J  ■' 

■u  I  _  -H '  •  'h*  -  h  ’  •  h'  =  b' 

i^j+1  ®  '  2^j+i  ®-  .,  2^j  '  3  J+1  .3  d  ■• 


Combining. these  terms  and  rearranging  them 


=  1  -  ^  (1  -  cos  kA) 

./  J  .1.  T  /o  1  ^ 


d+i/2. 


1/3.-  E^(|),-  |e3  (A)j  -V 


where 


oo 


n=l 


cos  knA  -  cos  k(n+l 


l)  A  I  E2(nA>.-  E2(n+1)  A 


The  terms  in  the  sumraation  can' be  reduced  to -a  single  integral  by  inter¬ 
changing  summation  and  integration. 

^  Hdji  |(1-  e  (e^^,  -  1) 

n=l  I  . 

+  (1  -  e"  -  1)  ^ 


',(26) 
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(I  -  e  M)  (.tto  ■  1)  ^ 


^ld|l 


;)  -  1 


-  e  “)  -  1) 

*  S>  -  1 


.1 


■yr 


cos  kA  (l  -  cos  kA)  sinh  --  +  sin^  ki^  (cosh  —  -  l) 
_  _ li  _ u 


■  cosh  —  -  cos  kA 


The  quantity  ^  is  a  function  of  the  two  dimensionless  quantities  kA  and  A 
in  which  kA  plays  the  role  of  a  parameter  to  be  varied  in  order  to  find  the 
least  stable  value. 


which 


First,  stability  is  considered  for  the  limiting  case’ of  A::?;^!  for 


y  =  -g  (l  f  2  cos  kA)  (l  -  cos  kA)  (a),  aV"^! 


In  this  case,  ^1  since  |J]|  is  less  than  the  first  term  insuring  that 

the  coefficient  of  ig  positive.  It  is  easy  to  see  that  the  most 

•  ^1/2  •  , 

stringent  stability  condition  is  realized  when  kA  =  ,(2n+l)n  ,  n  =  0,1,2  ... 


for  which 


where" 


f(A)  =  1/3  -  Ej^  (|)  -  I  (A)-  .  ■■■■.••• 

For  A':^-="l  ,  f(A)  =  1/3  and  ^  - — :>0  giving  the  limiting  stability  condition 


AfiL 


3  p-  Cy’.Ax  ^ 


8  ca  T 


3  ■ 


■1  .  . 


(27) 


As  is  expected  from  the  fact  that  the  difference  equation  becomes  just  the 
explicit  diffusion  equation  in  this  case.  Equation  (27)  is  the  diffusion 
stability  condition.  In  general,  f.(A)  is  smaller  than  1/3  as  shown  in 
Figure  1.  For  A  an  expansion  of  f  yields: 
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FIGURE  I. 
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where 


•  f  y  /n  (rA)  , 


/n  7  =  .5772 


Havinc  obtained  the  stability,  condition  for  Isirge  A  it  is  now  necessary 
to  investigate  smaller  values  of  A.  In  terms  of  defined  qua^ntities  the 
ajnplification  factor  is: 


oAt  , 
e  =  1  - 


'j+1/2 


—  (1  cos  IcA)  f  (a) 


-I,]  ■ 


While,  for  lare,e  A,  the  term  in'  ^  could  be  neglected  it  must  be  considered 
in  other  cases.  Some  of  the  properties  of  Y. 

1.  5]  (2jrn  +  kA)  =  ^  (hA)-,  n  =  1,2,3  ■... 

2.  Y;(2k  -  kA)  =  •  y;  (kA)  , 


y](kA  =  0)  =  0. 


5.  ^(kA  =  n)  =  4j 


A  •  ,  A 

-  •—  cosh  —  -  1 

e.  ^  -  tidu.il  0 

cosh  — 

-  —  sinh  — 

^  _ \L. 


cosh  —  +  1 
h 


lidu  7  0 


•  The  stability  behavior  for  A  .^1  is  more  complicated  mathematically 
than  the  diffusion  limit  discussed  above  because  Y  contains  the  leading 
term  in  Equation  (26).  It  is  also  true  that  Y  docs  not  take  a  simple  limit¬ 
ing  form  in  the  neighborhood  of  kA  =  0.  With  neither  a  closed  form  expression 
.for  ^]  in  the  general  case  nor  even  one  for  the  investigation  has 

taken  the  form  of  evalmting  integrals  which  bound  Y  from  above  and  below. 

For  A -1-^1  the  result  of  the  bounding  argument  indicates  that  the  coefficient 

of  — -  remains  positive,  that  it  has  a  maximum  value  near  kA  =  0,  and  that 

^i+1/2 

it  rapidly  falls  to  zero  at  kA  =  0.  The  stability  condition. 
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A 

’'j+1/2 


A  ^  -^1  , 


gives  the  limiting  time  interval 


-i 

2  <r  ca  T"^ 

This  value,  depending .not  at  .all  on  the  zone  size  Ax,  is  just  the  condition 
that  the  energy  removed  from  an  optically  thin  zone  in  the  time  interval 
not  drive  the  material  temperature  negative  enough  to  produce  instability. 

.The  numerical  .value  of  the  coefficient,  however,  depends  on  the  linearization 
of  the  equations .  A. value  representing  the  condition  tlmt  the  temperature 
in  the  actual  equations  must  not  become  negative  in  a  time  step  is: 


P.C^ 

At  £ - ^  •,  A  ^^1 

ct-  ca-  T;^ 


(28) 


which  is  to  be  preferred  over  the  previous  value . 

As  a  basis  for  stability  estimates  to  be  tested  with  sample  calculations 
the  two  limiting  values  given  above  are  incorporated  into  a  simple  inter¬ 
polation  formula : 


A- 

^  V 


(29) 


In  an  actual  calculation  in  which  the  quantities  entering  Equation  (29) 
depend  on  the  zone,  the  minimum  value  of  At  resulting  from  testing  all  zones 
of  the  problem  should  be  used  to  govern  the  time  step. 
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